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ABSTRACT 

The  problem  of  solving  the  equation  for  uniform  flow  in  an  open  channel 
is  reduced  to  that  of  evaluating  three  simple  integrals.  Practical  methods  of 
computing  these  integrals  are  discussed. 
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[2],  [3],  and  [4]  can  be  eubetituted  into  [1]  to  yield 


dx 


i  -  (yc/y)M 

So  -  {Q2/Ck)y‘k 


dy 


[5] 


where  the  critical  depth  yc  ie  given  by 


yc  =  (zc2/cz)1/M  .  [5'] 

In  moat  applications  [5]  is  integrated  to  give  the  distance  0  between  two 
depths  y^  and  y., : 


6  = 


1  -<ye/y>M 
So  -  (Q2/ck)y'N 


dy 


yz  *  yi 


[6] 


The  only  mathematical  consideration  to  be  observed  in  taking  the  integral 
of  [6]  is  that  no  zero  of  the  denominator  of  the  integrand  be  permitted  to  occur 
over  the  range  (y^  y2>*  .  [6]  can  be  given  the  more  convenient  form 

5  =  B(M,  N;  So,  yc,  Q2/Ck;  y2)  -  B(M,  N;  Sq,  yc,  Q2/Ck;  )  [7] 

where 

2  .  f  1  "  (yc/y)M  r  ! 

B(M,  N;  S  y  Q2/^;  y)  =  \  - f-  -  ■  ,N  dy  .  [7'] 

0  C  ^  J  SQ  -  (QZ/Ck)y  N 

It  is  convenient  to  study  the  function  B  . 


4i 

Physically, this  would  correspond  to  prohibiting  the  flow  profile  in  a  uniform 
channel  from  smoothly  crossing  the  normal  depth.  This  is  in  accord  with 
experimental  fact. 
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reduces  [9]  to 

Bu(y)  =  •  idr  J 


1  +  N  w  M  -  1  -  N 

k2  ^  "  (yc/yn  >  X2  " 

- rr t: - 


IT 


d\,  [13] 


Class  UI:  H2,  H3. 

In  this  class  S  =  0.  Thus  [7']  reduces  to 


oo  >  y  >  y  >0 
7  7n 


Bm(y)  =  — j— ^ -  f 

ni  (Q2/Ck)  J 


N  M 

y  -  yc  y 


N-M~J 


dy  0  <  y  <  oo 


Class  IV:  A2,  A3. 

In  this  class  Sq  <  0  and  [7*]  reduces  to 

i  r  1  -  (y^/y)1 


Biv  y)  =  S~  J  ~ 

o 

i  +  (ya/y) 

where  the  adverse  depth  y  is  given  by 

cl 

Q2 

ya 

5oCk 

Suppose  y  /y  ^  1  .  Define 

cl 

X3  " 

(ya/y)"N 

dy 


1/N 


Then 


(y./y) 


-N 


Biv(y)=  IT#"  I 


l/N 


-  (y,/yJM  K  “W 


1  -  M 


3  ''C'a. 

- ITT 


dX. 


[14] 


[15] 


[15-] 


[16] 


[17] 


0  <  y  ^  ya  <  oo 
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Suppose  ya/y  <  1  .  Define 


*4  •  <ya/v)N 


Then 


BIV(*>  =  BIV(>'a)  •  Tjl- 


r- 


M-l-N 


■4 


ya  ^  y  <  oo 

Claes  V:  HI,  Al;  C2. 

HI  and  Al  are  not  physically  realizable  (Ref.  1,  Ch.  9).  C2  is  the 

trivial  case  y  =  y  =  y„  ;  it  is  unstable, 
c  n 

The  most  interesting  fact  about  the  integrals  thus  far  presented  is 
that  they  can  be  simply  expressed  in  terms  of  three  fundamental  functions, 
the  backwater  integrals: 


•>>  ■ ) 

<8  y(z)(z)  =  J  zy  '  1 


dz  0  <  z  <  1 


0  <  z  <  oo 


[20a] 


[20b] 


■I 


One  can  thus  reduce  the  B's  to 


0  <  z  <  1 


BI(y)  =  "  Trk"  I  ®1+N  ®/-M+N  (al* 


[20c] 


al  =  (yn/y)' 
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BIlW  *  -  TO-  «  -V/N  <‘2>  •  (ve/v„)M  «  (m’.  1  <«2> 


*  <y„/y>  N 


l‘3»  -  VcM  «N-M+1  '“3» 


Q,  =  y 


Za  30(3) 


Biv(y)  =  iri-  ^I+n  (Q4)-(yc/ya)M  ®i(-m+n  (a4) 


a4  =  (ya/yrN 


ya  (^>-<ya/ya)M  ^  t24"] 

°  ”tt  ~ir 


In  the  next  section  methods  of  eva.lua.ting  the  (8'*'  (z)  will  be  discussed. 


3.  Evaluation  of  the  Backwater  Integrals 
The  evaluation  of  the  backwater  integrals  for  all  values  of  7  presents 
no  intrinsic  difficulty.  However,  it  can  become  rather  involved,  especially 
for  7  =  -m(m  =  0,  1,  2,  .  .  . ),  and  presents  practical  problems  in  computation. 
For  the  purposes  of  this  study  it  is  sufficient  ,  for  (B^  and  (B^,  to  consider 


only  those  portions  of  the  7-axis  covered  by 


1+N  1-M+N 


IT  > 


*  ;  or,  for  <8  ^  ,  by  1+N  and  1  -M  +  N  .  One  can  first  note  that 


and 
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The  restrictions  placed  upon  y  by  [26]  and  [27]  preclude  the 
existence  of  logarithmic  terms  in  any  of  the  <&^\z).  Hence,  it  is  convenient 
to  expand  the  integrands  of  [20a]  and  [20c]  in  infinite  series  and,  uniform 
continuity  being  seen  to  obtain,  to  integrate  these  series  term  by  term: 

«t‘>n  .  t28] 

n=o 
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It  happens  that  the  series  in  [28]  can  be  summed  at  once  (2)  to  yield 

{ 3 )  _Y 

(z)  =  —  F  T»  7+1;  ±«)  [29] 

where  F  is  an  ordinary  hypergeometric  function.  This  fact,  however,  is 
of  negligible  value  computationally  and  it  is  better  to  sum  [28]  by  less  direct 
methods. 

[28]  can  be  reduced  at  once  to 

( 1  ) 

<8y3  (z)=  zy  +Sy(  +  z)1  [30] 

where 
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The  Eq.  [32]  can  itself  be  simplified  by  noting  that 

y  41C—  =  z-y  f  ?  ii^L(z.)^-1dz 

JA  “ <n  +  *>  o  n=  1  n 

z 

=  z'7  y  (z')^"*  Lij^(+_z ')  dzf 
o 

Let  a  polynomial  approximation  to  Li^  be  formed: 

J 

UK(iz)=l  KCj(-zJ  +  KRJ-  (z) 
j=l 


Then 


O  £t 

s7(+z)  ^  oCr  J Ty~  zj  +Z'Y  J(*')ir;lRJ(±){»,)d»'  K  =  0  t35aJ 

j=l  o 

■I  (-T)k'1Lik(t*)+(-i')K-1X  Kcr-iV2i 

k=l  j=l 

z 

+  (*^)K*1  KRr  <z)+(-^-y  J (z')Y' W-  (z,) dz'  k>i.  t35bi 


The  choice  of  suitable  K  and  J  is  not  simple,  since  a  balance  must 
be  drawn  between  the  difficulties  of  computing  the  polylogarithm  functions 
as  opposed  to  those  of  summing  series.  In  addition,  S^(+z)  has  a 
singularity  at  z  =  1  which  adversely  affects  the  convergence  of  the  series 
which  obtain  for  K  ^  0  and  makes  it  difficult  to  obtain  suitable  approximating 
polynomials. 
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To  insure  that  the  integrals  in  [33]  were  well  behaved  it  was  decided 


Lik(+*) 


(+) 


to  approximate  — —  and  then  to  represent  j^Rj— 1 '  (z1)  by 


KRJL+)  <•'> 


*•  keP  <-■) 


[36] 


(+)  LiK(±*> 

where  j^Ej-'(z)  is  the  error  in  ^  +  .  If  then,  the  maximum  value  of 

(z)  |  over  (0,  z)  is  ^Mj^- (z)  the  error  terms  are,  in  absolute  value, 


less  than 


for  [35a]  and  less  than 


■frr  o1^2' 


I  |K-1 

I  y  I 

T+T- 


(z) 


[37a] 


[37b] 


for  [35b]  . 

The  construction  of  the  polynomial  approximations  is  discussed  in 
detail  in  Appendix  I;  it  is  sufficient  to  note  the  results  here.  For  integrals 
of  the  third  kind  K  =  1  was  found  to  be  adequate: 


w  (-) 

1M9 

< 

5  x  10 

-7 

[38.  0] 

9c  <-> 
11 

= 

-1.0000 

0000 

[38;  1] 

9.  (-) 

1  C2 

= 

0.4999 

8946 

[38.  2] 

9  c  (-) 

1  C  3 

= 

-0. 3330 

7718 

[38.  3] 

9  c  (-) 

1  C4 

= 

0. 2476 

3567 

[38.  4] 

9C  (-) 

1 C  5 

= 

-0. 1883 

5208 

[38.  5] 
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9c<-> 

1  c6 

=  0.  1315 

6190 

[38.6] 

9c<-> 

1  C7 

=  -0.0729 

0260 

[38.7] 

9c:(-> 

1  C8 

=  0.0265 

6504 

[38.8] 

9  c(-) 

1  C9 

=  -0.0045 

6742 

[38.9] 

Thua  , 


<B(t3>(z)  =  zy 


—  +Y  9  c  H  -X-  zj+  ^(SxlO^x-^-) 
y  L  i  j  j  +  7  y+T 

L  j=l 


For  the  integrals  of  the  second  kind 


[39] 


<8  (2»(z)  =  _L 
y  y 


[40] 


For  integrals  of  the  first  kind,  unfortunately,  it  proved  necessary  to  go  to 
K  =  4  to  obtain  a  polynomial  approximation  which  converged  with  sufficient 
rapidity: 


4M^+)  <  7xl0'7x(l -fn  (1  -  z)) 


[41.0] 


(This  odd  weighting  was  chosen  so  that  the  relative  accuracy  of  the 
approximation  would  not  be  greater  than  necessary  near  the  singularity  of 
Li^iz)  =  in(l-z)  at  z  s  1.  ) 


^  >0 

o 

II 

1.0000 

0046 

[41.1] 

9  c{+)  . 

4  c2 

0.0624 

1793 

[41.2] 

*>0 

O 

II 

0.0143 

7767 

[41.3] 

9  c(+)  . 

4  C4  ' 

-0.0150 

2213 

[41.4] 
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9  c(+)  _ 

4  C5  “ 

0.0873 

5126 

[41.5] 

9  c(+)  . 

4  C6  ' 

-0. 2087 

8801 

[41.6] 

9  c(+)  _ 

4  7  " 

0.  2831 

2643 

[41.7] 

9  c<+)  .. 

4  C8  ' 

-0.  1984 

0826 

[41.8] 

9  c<+>  - 
4  9  ' 

0.0572 

7062 

[41.9] 

Hence, 


<S  <1,(Z)  =  2^ 
y 


-i-  +  Lix  (z)  -  yLii^z)  +  yZ  Lij  (z) 


si  M+,Tt 

j=i 


+  Q'flx  10'7  x  x  z  X  (1  -  in(l-z))^ 


[42] 


Since  the  polylogarithms  have  been  extensively  tabulated(4),  the 

forms  [39],  [40],  and  [42]  provide  a  convenient  way  of  computing  the  ($^(z) 

by  hand.  The  attainable  accuracy  will,  of  course,  vary  with  the  choice  of  f  , 

y  ,  and  z:  it  is  believed  that  if  eight  figures  are  carried,  the  relative  error 

5 

in  the  answer  will  normally  be  less  than  1  part  in  10 

The  Eqs.  [39]  and  [40]  are  well  suited  to  use  with  digital  computers. 


[42]  is  not  quite  as  convenient  for  machine  use  since  it  requires  the  computation 
of  three  polylogarithm  functions  which  can  in  turn  necessitate  the  computation 
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o f  five  logarithms  plus  the  summation  of  two  slowly  convergent  series: 

(3) 

it  is  advisable  to  compute  <g  '  '  (z)  directly  from  [30]  when  the  series  for 
S^(z)  converges  with  suitable  rapidity.  From  [30'] 


oo 

n 

S  (z)  =  )  -4- 
y  / J  a  +  1 

n=  1 


I  *•  I  <  1 

-1/2  <7<  3/2 


[43] 


To  determine  the  z-range  over  which  [43]  is  useful  it  is  convenient  to 
rewrite  it  as 


sv(«> 


*  I  .-£? +  tn 

n=l 


(z) 


[44] 


where  the  truncation  error  T^(z)  is  given  by 


rN<*>  *  I  TTy 


n=N+l 

An  application  of  Cauchy's  integral  method  to[44']yields 


[44'] 


oo 


|tn(z)U  Mn  £ 

n=l 


&  *  UIN  f 


TN  -  l  )  +  t 


dt 


[45] 


which,  subject  to  the  substitutions 

|z|  =  e~p 
(N-l)  =  \ 


p  >  0 


[46a] 

[46b] 
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becomes 


Hence, 


where 


IV'I  ^*'NP  J  -T TT-  d* 


TN(z)|<e‘p  Ci(-Xp) 


Si(u) 


w 

•S  4- 


1*1  <i 

p  >  o 


Thus,  the  relative  error  (Slj^(z)  in  the  sum  is,  for  z  0  ,  bounded  by 

i<*N(*)|  £  “§■  I  &  (-^p)  I  •  [49] 

Taking  j(ft^(z)  |  <  5  x  10~°  (X.p  =  15), one  sees  that  N  =  64  will  suffice  for 

z  ^  0.75  .  Hence,  for  use  on  digital  machines  it  is  suggested  that  the  formula 

64  n  0  <  z  ^  3/4 

(B^tz)  =  zY  ^  nZ+T  -  -y  ^  7^  3/2  [50] 

n=°  7  /  0 

be  employed  when  possible. 

Rather  than  use  the  methods  outlined  above  to  produce  extensive 
tables  of  the  <B,^(z),  which  would,  in  this  era  of  digital  machinery,  be 
obsolescent  from  birth,  the  equations  [39],  [40],  [42],  and  [50]  were 
combined  to  produce  a  single  subprogram  for  the  IBM  7090  computer. 

Copies  of  this  subprogram  can  be  obtained  by  writing  the  author;  a  complete 
description  of  it  is  given  in  Appendix  II. 


i 
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4.  Conclusions 

It  was  shown  that  the  problem  of  integrating  the  differential  equation 
governing  the  uniform  flow  of  liquid  in  an  open  channel  can  be  reduced  to  the 
evaluation  of  three  transcendental  functions.  This  latter  problem  was  shown  to 
reduce  to  that  of  computing  certain  simple  algebraic  quantities  and  calculating 
various  low-order  polylogarithm  functions. 

No  tables  of  backwater  integrals  were  presented,  since  it  was  concluded 
that  -  in  view  of  the  extensive  and  rapidly  growing  use  of  digital  computers  - 
such  tables  would  not  be  of  widespread  interest.  Instead,  a  subprogram 
was  written  for  use  on  an  IBM  7090.  With  the  aid  of  this  subprogram  the 
engineer  can  rapidly  construct  a  program  to  calculate  whatever  flow  profile 
he  desires. 

An  area  where  considerable  future  progress  can  be  made  is  that  of  the 
derivation  of  approximating  polynomials  for  the  several  polylogarithm  functions. 
A  not  unrealizable  goal  should  be  the  polynomial  approximation  of  Li^ (z) 
which  would  reduce  the  evaluation  of  <8^  (z)  to  the  summation  of  a  rather 
long  series  and  the  computation  of  the  logarithm. 
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APPENDIX  I 


Construction  of  the  Polynomial  Approximations 
The  problem  to  be  dealt  with  here  is  that  of  properly  approximating 
certain  transcendental  functions.  It  was  solved  as  follows. 

Let  f  ( x )  be  a  function  which  is  to  be  approximated  by  a  polynomial 


Pj(x)  over  an  interval  (0,  1) .  Let 


Dj(x)  =  f (x)  -  Pj(x) 


[A-l] 


Further,  let  g(x)  be  some  positive  definite  function  over  (0,  1),  and  let 


Zj(  x )  =  Dj(x)/g(x) 


[A -2] 


be  the  error  relative  to  the  normalizing  function  g(x).  The  problem  is  to 

adjust  the  coefficients  a.  of 

J 


Pj(x)  -  2  .jX-i 
j=0 


[A-3] 


so  that  Zj(x)  falls  within  some  desired  bounds. 

The  method  used  is  based  upon  that  of  Spitzbart  and  Macon  (4)  .  Let 
Zj(x)  be  specified  over  a  set^x^j.  of  J  +  2  distinct  points  as 


ZT(x  )  =  \  d 
J  p  p 


[A-4] 


where  the  have  been  specified  and  d  is  some  constant  to  be  determined. 
Then 

f(x  )  -  PT(x  )  =  d\  g(x  ) 

'  p  '  J  p  PBP 


[A -5] 
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d  ,  a  measure  of  the  goodness  of  fit,  and  the  coefficients  a.  are  determined 


J 


by  solving  the  system  [A -5]  .  Let 

/:•  \ 

A  = 

by 


F  = 


W  = 


r 


\ 


£<xo)  ^ 

£(xJ+l^ 

1  X  • 

o 

J 

•  •  X 

0 

g(Xo> 

1  x‘  ' 

J 

•  •  X1 

g(*! ) 

S 

1  XJM 

’  J 

XJ+1 

8(xJ+l) 

XJ+1  / 

[A -6a] 


[A -6b] 


[A-6c] 


With  this  abbreviated  notation  [A- 5]  can  be  rewritten  as 


WA  =  F  [A-7] 

[A-7]  will  have  a  non-trivial  solution,  i.  e.  ,  W"1  will  exist,  if  and  only  if 
no  polynomial  of  degree  less  than  J+l  passes  through  the  points 
( xp,  ^  g(  xp) )  (p  =  0,  1,  2,  .  .  . ,  J  +  1 )  . 

For  the  approximations  desired  here  this  condition  was  met  by  setting 


X.  = 
P 


(-1)P 


(p  =  0,  1,  .  .  . ,  J) 


[A-8a] 
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and 

XJ+1  =  XJ  • 


[A-8b] 


Of  course,  the  formalism  leading  to  [A-7]  does  not  really  solve  the 
approximation  problem  completely, since  it  does  not  specify  either  g(x) 
or  the  set  |xpj  •  The  choice  of  g(x)  was  left  arbitrary  and  was  varied  from 
one  f(x)  to  another.  The  ^x^j  problem  was  solved  by  requiring  the  extremals 
of  Zj(x)  to  coincide  with^x^  .  For  monotone  f  (x)  this  leads  to  the 
danger  of  having  only  J  +  1  extremals,  and  in  this  case  Xj+j  was  chosen  to 
be  1,  Xj  being  unchanged  if  x  =  1  wad  an  extremal,  or  being  placed  between 
x  =  Xj  j  and  x  =  1  if  x  =  1  was  an  ext  remal.  This  choice  of£x^  meant  in 
practice  that  the  set  £x^j  was  found  iteratively:  a  set  jx^^ j  was  chosen,  the 
Sj  found,  and  a  curve  of  Z^(x)  computed;  from  this  Z^  a  new  set  of 
extremals,  i.e.,  was  found  and  the  process  repeated;  the  iteration 

was  continued  in  this  manner  until  convergence  resulted  or  until  Z^(x) 
and/or  d  became  so  small  that  computing  error  rendered  further  iteration 
unprofitable.  This  method  of  selecting  the  and  the  ^x^j  effect  adjusted 
D^(x)  so  that  its  extremals  were  of  alternating  sign  and  so  that  |g(x)d  | 
was  an  envelop  for  it. 

A  FORTRAN  program  embodying  these  features  was  written  and  run 
on  the  IBM  7090  computer  at  the  Harvard  University  Computing  Center.  It 
was  found  that  the  W-matrices  tended  to  be  fairly  poorly  conditioned  with  the 
result  that  it  was  difficult  to  procure  solutions  for  [A-7]  when  J  exceeded  ten. 
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Since  <B  ^‘{z)  is  singular  at  z  =  1,  the  weighting  function 
g  (x)  =  1  -  in  (1-z) 


[A-9] 


was  chosen,  and  approximations  to  Li^(z)/z  for  various  J  and  p  =  2,  3,  4  were 

generated.  The  results  of  these  calculations  are  summarized  in  Fig.  I  . 

th 

The  singularity  at  z  =  1  in  the  p  derivative  of  Li^(z)  is  seen  to  affect 

powerfully  the  goodness  of  fit.  The  case  p  =  4,  J  =  8  was  chosen  as  adequate 

and  the  results  presented  in  [41]  were  adopted.  With  an  approximating  program 

of  higher  accuracy  it  seems  likely  that  sufficiently  accurate  results  could  be 

obtained  with  p  =  3,  J  >  16  .  It  seems  probable  that  the  obtaining  of  a  good 

approximation  for  p  =  2  will  require  a  quite  large  value  of  J. 

(3) 

The  fitting  of  (B^  '  (z)  was  much  simpler, since  it  and  its  derivatives 
possess  no  singularity  over  (0,  1).  In  this  case  a  weight  function  equal 
approximately  to  unity  was  chosen  and  the  results  quoted  in  [38]  obtained 
for  p  =  1,  J  =  8  . 
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APPENDIX  II 

A  Program  lor  Computing  the  Backwater  Integrals 

A  FAP  coded  FORTRAN  type  subprogram,  named  BKW,  was  written 
to  compute  the  backwater  integrals.  It  was  subsequently  debugged  and  tested 
on  the  IBM  7090  at  the  Harvard  University  Computing  Center.  Binary  decks 
for  this  routine  may  be  obtained  by  writing  the  author;  these  decks  also 
contain  the  special  subroutines  AGLOG,  AGPWR,  and  AGTLiG  called  by  BKW. 

A  listing  of  BKW  is  presented  in  Figs.  Ila  through  lie.  The  program 
for  computing  <B  ^(z)  is  seen  to  be  divided  roughly  into  five  major  sub¬ 
sections.  In  the  first  subsection  the  program  decides  which  of  the  three 
integrals  (I  =  1,  2,  3)  is  to  be  evaluated  and  transfers  control  to  the  appropriate 
subsection.  There  is  a  separate  subsection  (Bl,  B2,  or  B3)  for  each  value 
of  f  .  Finally, there  is  a  subsection  (ERR  and/or  RET)  which  controls  returns 
to  the  calling  program. 

Bl.  This  subsection  first  examines  z  :  for  z  <  0  or  z  ^  1,  an  error 
return  (transfer  to  ERR)  is  specified,  for  0  <  z  ^  0.  75,  an  evaluation  by  the 
series  [50]  (transfer  to  SERI)  is  specified)  for  0.75<z<1.00,  an 
evaluation  by  means  of  E<j.  [42]  is  specified.  The  first  step  in  the  evaluation 
by  either  [50]  or  [42]  is  the  examination  of  y  :  if  the  y  specified  equals  a 
value  held  in  storage  (zero  if  this  is  the  first  pass  through  SERI  or  the 
corresponding  routine  for  large  z  ,  otherwise  the  value  of  y  on  the  last 
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pass  through  the  appropriate  routine),  the  values  of  (or  ^  cj+^  - ) 

are  known  to  be  computed  and  held  in  storage  already  so  that  <8^M(Z)  can 
be  computed  at  once  from  [42]  (or  [50]  )  and  control  transferred  to  RET;  if  y 
does  not  agree  with  the  value  held  in  storage,  it  is  checked  to  make  sure  that 
it  is  within  the  prescribed  range  of  values  and,  assuming  that  it  passes  this 
test  (transfer  to  ERR  otherwise),  the  appropriate  constants  are  calculated  and 
the  integral  evaluated  as  before.  In  the  course  of  these  calculations  various 
special  subroutines  are  utilized:  AGLOG,  a  27  bit  accuracy,  floating  point 
logarithm  routine;  AGPWR,  a  26  bit  accuracy,  floating  point  power  routine; 
AGDJLG,  a  routine  for  Li^fz);  AGTLG,  a  routine  for  Li^  (z). 

B2.  This  subsection  first  obtains  y  and  z  and  tests  to  see  that  they 
lie  within  the  allowed  ranges.  Failure  of  thiB  test  causes  program  control  to 
be  transferred  to  ERR.  Success  of  this  test  results  in  the  evaluation  of 
<8^(z)  directly  from  [40];  control  is  then  transferred  to  RET. 

B3.  This  subsection  is  roughly  similar  to  that  for  (8  ^^(z),  except 
that  the  integrals  for  all  allowable  z  are  computed  from  the  same  formula. 
First,  z  is  checked  for  range.  Then  y  is  checked,  and,  if  necessary,  the 
constants  J  cj  ^  are  recomputed.  Finally,  <B^  (z)  is  computed  from 

[39],  and  control  is  transferred  to  RET. 
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ERR  and  RET.  ERR  is  an  error  routine  which  clears  the  accumulator 
and  turns  on  sense  light  4;  it  then  transfers  control  to  RET.  RET  resets  the 
three  index  registers  and  returns  control  to  the  calling  program. 

*  *  * 


USAGE.  The  appropriate  calling  sequence  in  FAP  is 
A  TSX  BKW,  4 

TSX  L 

TSX  GAMMA 


TSX 


(i), 


The  normal  return  is  to  A+4  within  (8^  '(z)  in  the  accumulator;  in  case  of  an 
error  in  the  range  of  z  or  y  ,  return  is  to  A+4  but  the  accumulator  will  be 
clear  and  sense  light  4  will  be  on.  (Note  that  L  must  be  a  FORTRAN  integer 
and  that  GAMMA  and  Z  must  be  FORTRAN  floating  point  numbers.  )  The 
appropriate  calling  sequence  in  FORTRAN  is 

Y  =  BKW ( L,  GAMMA,  Z) 
where  Y  is  the  value  of  the  desired  integral. 

TIMING.  This  is  quite  variable.  It  runs  from  about  10  ms  for  the 

worst  cases  of  <8^(z)  down  to  2+  ms  for  <8^^(z)  and  2  ms  for  (8^(a). 

5 

ACCURACY.  Error  is  expected  to  be  less  than  1  part  in  10  in  most 
cases.  For  <B  (0.  75  <  z  <  1. 00)  it  is  normally  less  than  1  part  in  10^. 


ENTRY 

BKW 

SXO 

#-2.4 

SXD 

ERAS+1 t 1 

SAVE  INDICES 

SXO 

ERAS+2.2 

SXD 

ERAS+4.4 

CLA# 

1.4 

GET  L 

CAS 

CNST 

TRA 

*+2 

TRA 

B1 

L  =  1 

CAS 

CNST+1 

TRA 

#+2 

TRA 

B2 

L  =  2 

CAS 

CNST  +2 

TRA 

*+2 

TRA 

B3 

L  =  3 

TRA 

ERR 

CLA* 

2.4 

GET  AMD  SAVE  GAMMA 

STO 

ERAS 

CLA# 

3,4 

GET  Z 

TMI 

ERR 

Z  NEGATIVE 

T2E 

ERR 

Z  ZERO 

CAS 

CNST+3 

TRA 

ERR 

Z  TOO  LARGE 

TTR 

*-l 

STO 

Z 

Z  OK.  SAVE 

CAS 

CNST +4 

TRA 

#+6 

Z  GT  0.75 

TTR 

#+l 

TSX 

SERI. 4 

HTR 

ERAS 

HTR 

Z 

TRA 

RET 

CLA 

ERAS 

CAS 

G1 

TTR 

#+2 

TRA 

#+2 

OLD  CONSTANTS  SUFFICE 

TRA 

CONI 

GET  NEW  CONSTANTS  FOR  ASYMP  EVAU 

CLA 

CNST+3 

BEGIN  ASYMPTOTIC  EVALUATION 

FDP 

G1 

STO 

TOT 

CLA 

CNST  +  3 

FSA 

Z 

STO 

ERAS 

TSX 

5AGL0G.4 

HTR 

FRAS 

HTR 

ERAS 

CHS 

FAD 

TOT 

STO 

TOT 

LI1  INCLUDED 

TSX 

SAGDLG.4 

HTR 

Z 

XCA 

FMP 

G1 

CHS 

FAO 

TOT 

STO 

TOT 

L  I  2  INCLUDED 

TSX 

SAGTLG.4 

HTR 

Z 

FI6.no  LISTING  OF  BKW  (PACE  I  OF  5) 


XCA 

FMP 

XCA 

FMP 

FAD 

STO 

AXT 

PXD 

FAD 

XCA 

FMP 

TXI 

TXL 

XCA 

FMP 

XCA 

FMP 

XCA 

FMP 

CHS 

FAD 

STO 

TSX 

HTR 

HTR 

XCA 

FMP 

TRA 

SERI  SXO 
CLA* 
STO 
CLA* 
CAS 
TTR 
TRA 
TRA 
AXT 
PXD 
FAD 
XCA 
FMP 
TXI 
TXL 
FAD 
STO 
TSX 
HTR 
HTR 
XCA 
FMP 
LXD 
TRA 

CONI  TZE 
CAS 
TRA 
TTR 
TRA 
CAS 


G1 

G1 

TOT 

TOT 

0.1 

0.0 

Sl+8.1 

Z 

*+1.1.1 

*-4.1. 8 

G1 

G1 

Gl 

TOT 

TOT 

SAGPWR.4 

Z 

Gl 

TOT 

RET 

ERA3+3.4 

2.4 
Z 

1  .4 
G1S 
*+2 
*  +  2 
CON1S 
0.4 
0.0 

S1S+64.4 

z 

*+1.4.1 

*-4.4.63 

SIS 

ERAS 

SAGPWR.4 

Z 

Gl  S 
ERAS 

ERAS+3  .4 

3.4 
ERR 

CNST+5 

*+3 

»-l 

ERR 

CNST+6 


L I  3  INCLUDED 

START  APPROXIMATING  SERIES 


COMPLETE  SERIES 


B  1  COMPLETE 

SAVE  XR ( 4 ) 

SAVE  7 
GET  GAMMA 


OLD  CONSTANTS  SUFFICE 
NEW  CONSTANTS  NECESSARY 
BEGIN  SERIES  SUMMATION 


CHECK  RANGE  OF  GAMMA 


FlC.Hb  LISTING  OF  BKW  (PACE  2  OF  5) 


TRA 

ERR 

TTR 

*+l 

STO 

G1 

AXT 

0.1 

CLA 

CNST+7 

STO 

ERAS+5 

CLA 

ERAS+5 

FAD 

G1 

STO 

ERAS+6 

CLA 

ERAS+5 

FDP 

ERAS+6 

FMP 

C+8.1 

STO 

Sl+8.1 

CLA 

ERAS+5 

FSB 

CNST+3 

STO 

ERAS+5 

TXI 

*+l . 1.1 

TXL 

*-11.1.8 

TRA 

Bl+21 

CON  IS 

TZE 

ERR 

CAS 

CNST+5 

TRA 

*+3 

TTR 

*-l 

TRA 

ERR 

CAS 

CNST+6 

TRA 

ERR 

TTR 

*+l 

STO 

G1S 

AXT 

0.4 

CLA 

CNST+6 

STO 

ERAS+5 

CLA 

ERAS+5 

FAD 

G1S 

STO 

ERAS+6 

CLA 

CNST+3 

FDP 

ERAS+6 

STO 

SlS+64,4 

CLA 

ERAS+5 

FS3 

CNST+3 

STO 

ERAS+5 

TXI 

*+1.4.1 

TXL 

*-10.4,64 

TRA 

SERI +8 

ERR 

PSE 

100 

PXD 

0,0 

TRA 

RET 

RET 

LXD 

ERAS+1,1 

LXD 

ERAS+2 ,2 

LXD 

ERAS+4,4 

TRA 

4.4 

32 

CLA¬ 

2.4 

ST  0 

C-2 

t.;i 

ERR 

TZE 

ERR 

CAS 

C.MST+9 

TRA 

ERR 

TTR 

*+] 

CLA* 

1  ,4 

FIS.  He 


SAVE  NEW  VALUE  OF  GAMMA 
SET  UP  NEW  CONSTANTS 


RETURN  TO  B1  EVALUATION 
CHECK  RANGE  OF  GAMMA 


SAVE  NEW  VALUE  OF  GAMMA 
SET  UP  NEW  CONSTANTS 


RETURN  TO  SERI 
TURN  ON  SENSE  LIGHT  4 
CLEAR  AC 
PREPARE  TO  EXIT 


GO  SACK  TO  MAIN  PROGRAM 
OBTAIN  GAMMA  A  N  u  CHECK  RAN 


03T  A I Z  and  CHECK  RANGE 

LISTING  OF  BKW  (PAGE  3  OF  5) 


TV! 

ERR 

TZE 

ERR 

STO 

Z 

TSX 

SAGPWR.4 

Z**GAMMA 

HTR 

Z 

HTR 

G2 

FOP 

G2 

XCA 

TRA 

RET 

CLA# 

2«4 

GET  GAMMA 

CAS 

G3 

TRA 

*+2 

TRA 

*+2 

OLD  CONSTANTS  SUFFICE 

TRA 

COM3 

NEW  CONSTANTS  NEEDED 

CLA* 

3.4 

GET  AND  TEST  Z 

TMI 

ERR 

TZE 

ERR 

CAS 

CNST+3 

TRA 

ERR 

TTR 

#  +  l 

STO 

Z 

AXT 

0.1 

SUM  srpicc 

PXO 

0.0 

FAD 

S3+8 .1 

XCA 

FMP 

z 

TXI 

*+1.1.1 

TXL 

*-4* 1.8 

STO 

ERAS 

CLA 

CNST+3 

FOP 

G3 

XCA 

FAD 

ERAS 

STO 

ERAS 

FIND  93 

TSX 

SAGPWR.4 

HTR 

Z 

HTR 

G3 

XCA 

F<  ;p 

ERAS 

TRA 

RET 

TZE 

ERR 

TEST  GAMMA 

CAS 

CNST+5 

TRA 

*+3 

TTR 

*-l 

TRA 

ERR 

CAS 

CNST+6 

TRA 

ERR 

TTR 

»+l 

STO 

G3 

GAMMA  OK 

AXT 

0.2 

SET  UP  SERIES  REVISION 

CLA 

CM  ST. +7 

STO 

ERAS+6 

CLA 

ERAS+6 

REVISE 

FAD 

G3 

STO 

ERAS 

CLA 

D  +  8.2 

FDP 

ERAS 

fmp 

ERAS+6 

FIS.  Ed  LISTING  OF  BKW  (PAGE  4  OF  5) 


CLA 

ERAS+6 

FSB 

CNST+3 

STO 

ERAS+6 

TX  [ 

*+1.2*1 

TXL 

*-11*2.8 

T3\ 

33  +  5 

CNST 

PZC 

3*0*1 

PZE 

0.0.2 

PZE 

0,0.3 

DEC 

1.00 

DEC 

0.75 

DEC 

-0.50 

DEC 

1.50 

DEC 

9.00 

DEC 

64.00 

DEC 

7.00 

C 

DEC 

+0.100000046E+01 

L  1 4 

CONSTANTS 

DEC 

+0.624179326E-01 

DEC 

+0.143776655E-01 

DEC 

-G.150221251E-01 

DEC 

+0.873512559E-01 

DEC 

-0.208788007E-00 

DEC 

+0.283126429E-00 

DEC 

-0. 198408261 E-00 

DEC 

+0. 572706185E-01 

0 

DEC 

-0.999999993E+00 

LI1 

CONSTANTS 

DEC 

+0.499989465E-00 

DEC 

-0.333077 177 E-00 

DEC 

+0.247 63567 OE -00 

DEC 

-0.186352078E-00 

DEC 

+0.131561905E-00 

DEC 

-U.729O25990E-01 

DEC 

+0.265650421E-01 

DEC 

-0.456741825E-02 

ERAS 

BSS 

7 

Z 

BSS 

1 

61 

BSS 

1 

G  IS 

BSS 

1 

62 

BSS 

1 

63 

BSS 

1 

TOT 

BSS 

1 

SI 

BSS 

9 

SIS 

BSS 

65 

S3 

BSS 

9 

END 

FIG.  He  LISTING  OF  BKW  (PAGE  5  OF  5) 
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(3)  7 

For  it  is  normally  less  than  2  parts  in  10  and  often  less  than  2  parts 

8  (11 

in  10  For  <BV  (z)  (0.  00  <  z  ^  0.  75)  the  error  is  normally  confined  to  the 

^  (i\  g 

eighth  figure.  The  error  in  (8'  '(z)  is  normally  less  than  2  parts  in  10 

7  (1) 

Regions  of  high  relative  error  can  be  expected  in  (8^'(z)  for  y  <  0  and 
values  of  z  for  which  <8^(z)  =  0  ;  however,  even  here  several  good  decimal 
places  will  normally  be  available. 
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